Above Ground Biomass (AGB) Estimation

Loading Required Libraries

library(sf)
library(knitr)
library(BIOMASS)
library(ggplot2)
library(dplyr)

#old data
data_forAGB <- data_clean
#nrow(data_forAGB)
data_forAGB

#new data
new_data_for_AGB <- new_data_clean
new_data_for_AGB

#Checking and retrieving tree taxonomy
#old data

Taxo_old <- correctTaxo(
  genus = data_forAGB$Genus.species,
  useCache = TRUE, verbose = FALSE)

saveRDS(Taxo_old, file = "~/Remote-sensing/taxo_old_data.rsd")

data_forAGB$GenusCorrected <- Taxo$genusCorrected
data_forAGB$SpeciesCorrected <- Taxo$speciesCorrected

Taxo_new <- correctTaxo(
  genus = new_data_for_AGB$Genus_species,
  useCache = TRUE, verbose = FALSE)

saveRDS(Taxo_new, file = "~/Remote-sensing/taxo_new_data.rsd")

new_data_for_AGB$GenusCorrected <- Taxo_new$genusCorrected
new_data_for_AGB$SpeciesCorrected <- Taxo_new$speciesCorrected

#Retrieving Family and Order Information #APG (Angiosperm Phylogeny Group) family and order (optional)

APG_old <- getTaxonomy(data_forAGB$GenusCorrected, findOrder = TRUE)
data_forAGB$familyAPG <- APG_old$family
data_forAGB$orderAPG <- APG_old$order

APG_new <- getTaxonomy(new_data_for_AGB$GenusCorrected, findOrder = TRUE)
new_data_for_AGB$familyAPG <- APG_new$family
new_data_for_AGB$orderAPG <- APG_new$order

#Estimating Wood Density

wood_densities_old <- getWoodDensity(
  genus = data_forAGB$GenusCorrected,
  species = data_forAGB$SpeciesCorrected,
  stand = data_forAGB$Plot_name
)

# The reference dataset contains 16467 wood density values??

wood_densities_old

wood_densities_new <- getWoodDensity(
  genus = new_data_for_AGB$GenusCorrected,
  species = new_data_for_AGB$SpeciesCorrected,
  stand = new_data_for_AGB$Plot_name
)

wood_densities_new

#ดึงค่า meanWD จาก df wood_densities ไปใส่ใน df คอลัมน์ใหม่
data_forAGB$WD <- wood_densities_old$meanWD
new_data_for_AGB$WD <- wood_densities_new$meanWD

data_forAGB

new_data_for_AGB

# dataset ของเรา มีต้นไม้กี่ต้นที่ได้ค่า wood density จากแต่ละระดับ (species/genus/family)
sum(wood_densities_old$levelWD == "species")

[1] 7264

sum(wood_densities_old$levelWD == "genus")

[1] 195

sum(!wood_densities_old$levelWD %in% c("genus", "species"))

[1] 1072

# dataset ของเรา มีต้นไม้กี่ต้นที่ได้ค่า wood density จากแต่ละระดับ (species/genus/family)
sum(wood_densities_new$levelWD == "species")

[1] 7294

sum(wood_densities_new$levelWD == "genus")

[1] 195

sum(!wood_densities_new$levelWD %in% c("genus", "species"))

[1] 1067

#Checking Missing Values After Wood Density Assignment

#check missing value แต่ละคอลัมน์
colSums(is.na(data_forAGB))

       ExactDate        Plot_name          Quadrat              Tag          StemTag    Genus.species               QX 
               0                0                0                0                0                0                0 
              QY               PX               PY              DBH              HOM             Code   GenusCorrected 
               0                0                0                0               56             6298                0 
SpeciesCorrected        familyAPG         orderAPG               WD 
               0                9                9                0 

#check missing value แต่ละคอลัมน์
colSums(is.na(new_data_for_AGB))

       ExactDate        Plot_name          Quadrat              Tag          StemTag    Genus_species               QX 
               0                0                0                0                0                0                0 
              QY               PX               PY              DBH              HOM             Code       Note_Genus 
               0                0                0                0               56             1913                0 
  GenusCorrected SpeciesCorrected        familyAPG         orderAPG               WD 
               0                0                0                0                0 

#ดูว่ามีแถวไหนบ้างที่ missing
data_forAGB %>%
  filter(is.na(familyAPG) | is.na(orderAPG))

utm_conner <- read.csv("~/Remote-sensing/DATA/utm_corner.csv")
head(utm_conner)

conner_filtered <- utm_conner %>%
  filter(Plot_name == unique(data_forAGB$Plot_name) & X_loc == 0 & Y_loc == 0)

head(conner_filtered)

tree_data_old <- data_forAGB %>%
  left_join(conner_filtered, by = "Plot_name") %>%
  mutate(
    abs_x = X_UTM + PX,
    abs_y = Y_UTM + PY  
  ) %>%
  dplyr::select(abs_x, abs_y) 

tree_data_new <- new_data_for_AGB %>%
  left_join(conner_filtered, by = "Plot_name") %>%
  mutate(
    abs_x = X_UTM + PX,
    abs_y = Y_UTM + PY  
  ) %>%
  dplyr::select(abs_x, abs_y) 

tree_data_old

tree_data_new

old_tree_sf <- st_as_sf(tree_data_old, coords = c("abs_x", "abs_y"), crs = 32647)
new_tree_sf <- st_as_sf(tree_data_new, coords = c("abs_x", "abs_y"), crs = 32647)

old_tree_longlat <- st_transform(old_tree_sf, crs = 4326)
new_tree_longlat <- st_transform(new_tree_sf, crs = 4326)

coords_longlat_old <- st_coordinates(old_tree_longlat)
tree_data_old$lon <- coords_longlat_old[, 1]
tree_data_old$lat <- coords_longlat_old[, 2]

coords_longlat_new <- st_coordinates(new_tree_longlat)
tree_data_new$lon <- coords_longlat_new[, 1]
tree_data_new$lat <- coords_longlat_new[, 2]

coords_tree_old <- apply(tree_data_old[c("lon", "lat")], 2, mean)
coords_tree_new <- apply(tree_data_new[c("lon", "lat")], 2, mean)

Estimating Tree Height (H)

H_chave_old <- retrieveH(
  D = data_forAGB$DBH,
  coord = coords_tree_old
)

H_chave_new <- retrieveH(
  D = new_data_for_AGB$DBH,
  coord = coords_tree_new
)

data_forAGB$H_chave <- H_chave_old$H
new_data_for_AGB$H_chave <- H_chave_new$H

old_HD_model <- modelHD(
  D = data_forAGB$DBH,
  H = data_forAGB$H_chave,
  method = "log2",
  useWeight = TRUE
)

new_HD_model <- modelHD(
  D = new_data_for_AGB$DBH,
  H = new_data_for_AGB$H_chave,
  method = "log2",
  useWeight = TRUE
)

Computing Above Ground Biomass (AGB)

data_forAGB$AGB <- computeAGB(
  D = data_forAGB$DBH,
  WD = data_forAGB$WD,
  H = data_forAGB$H_chave 
)

new_data_for_AGB$AGB <- computeAGB(
  D = new_data_for_AGB$DBH,
  WD = new_data_for_AGB$WD,
  H = new_data_for_AGB$H_chave 
)

Propagating Errors in AGB Estimation (Monte Carlo)

D_error_prop_old <- AGBmonteCarlo(
  D = data_forAGB$DBH, WD = data_forAGB$WD, H = data_forAGB$H_chave,
  Dpropag = "chave2004",
  errWD = rep(0,nrow(data_forAGB)), errH = 0 
)

D_error_prop_new <- AGBmonteCarlo(
  D = new_data_for_AGB$DBH, WD = new_data_for_AGB$WD, H = new_data_for_AGB$H_chave,
  Dpropag = "chave2004",
  errWD = rep(0,nrow(new_data_for_AGB)), errH = 0 
)

new_data_for_AGB

WD_error_prop_old <- AGBmonteCarlo(
  D = data_forAGB$DBH, WD = data_forAGB$WD, H = data_forAGB$H_chave,
  errWD = wood_densities_old$sdWD,
  Dpropag = 0 , errH = 0 
)

WD_error_prop_new <- AGBmonteCarlo(
  D = new_data_for_AGB$DBH, WD = new_data_for_AGB$WD, H = new_data_for_AGB$H_chave,
  errWD = wood_densities_new$sdWD,
  Dpropag = 0 , errH = 0 
)

H_model_error_prop_old <- AGBmonteCarlo(
  D = data_forAGB$DBH, WD = data_forAGB$WD, 
  HDmodel = HD_model, 
  Dpropag = 0 , errWD = rep(0,nrow(data_forAGB)) 
)

H_model_error_prop_new <- AGBmonteCarlo(
  D = new_data_for_AGB$DBH, WD = new_data_for_AGB$WD, 
  HDmodel = HD_model, 
  Dpropag = 0 , errWD = rep(0,nrow(new_data_for_AGB)) 
)

error_prop_old <- AGBmonteCarlo(
  D = data_forAGB$DBH,
  WD =  data_forAGB$WD,        
  H =  data_forAGB$H_chave, 
  Dpropag = "chave2004",           
  errWD = wood_densities_old$sdWD,    
  errH = H_model_error_prop_old$sdAGB
)

error_prop_new <- AGBmonteCarlo(
  D = new_data_for_AGB$DBH,
  WD =  new_data_for_AGB$WD,        
  H =  new_data_for_AGB$H_chave, 
  Dpropag = "chave2004",           
  errWD = wood_densities_new$sdWD,    
  errH = H_model_error_prop_new$sdAGB
)

error_prop_old[(1:4)]

$meanAGB
[1] 2379.039

$medAGB
[1] 2378.968

$sdAGB
[1] 39.31173

$credibilityAGB
    2.5%    97.5% 
2303.432 2453.130 

sum(data_forAGB$AGB)

[1] 1839.656

error_prop_new[(1:4)]

$meanAGB
[1] 2405.423

$medAGB
[1] 2405.962

$sdAGB
[1] 41.77771

$credibilityAGB
    2.5%    97.5% 
2326.113 2487.241 

sum(new_data_for_AGB$AGB)

[1] 1845.819

old_AGB_by_plot <- summaryByPlot(AGB_val = error_prop_old$AGB_simu, plot = data_forAGB$Plot_name, drawPlot = TRUE)

new_AGB_by_plot <- summaryByPlot(AGB_val = error_prop_new$AGB_simu, plot = new_data_for_AGB$Plot_name, drawPlot = TRUE)

Plot Spatial Configuration and Subplot Division

Ranong_plot_corner <- utm_conner[utm_conner$Plot_name == "La-un (Ranong)",]

Ranong_plot_corner

fieldplot <- check_plot_coord(
  corner_data =  Ranong_plot_corner,
  proj_coord = c("X_UTM", "Y_UTM"),  
  rel_coord = c("X_loc", "Y_loc"),
  plot_ID = "Plot_name",
  trust_GPS_corners = FALSE,
  draw_plot = TRUE, rm_outliers = FALSE)

fieldplot_old <- check_plot_coord(
  corner_data =  Ranong_plot_corner,
  proj_coord = c("X_UTM", "Y_UTM"),  
  rel_coord = c("X_loc", "Y_loc"),
  plot_ID = "Plot_name",                      
  trust_GPS_corners = FALSE,
  tree_data = data_forAGB,
  tree_coords = c("PX", "PY"),
  tree_plot_ID = "Plot_name",           
  draw_plot = TRUE, 
  rm_outliers = FALSE
)

fieldplot_new <- check_plot_coord(
  corner_data =  Ranong_plot_corner,
  proj_coord = c("X_UTM", "Y_UTM"),  
  rel_coord = c("X_loc", "Y_loc"),
  plot_ID = "Plot_name",                      
  trust_GPS_corners = FALSE,
  tree_data = new_data_for_AGB,
  tree_coords = c("PX", "PY"),
  tree_plot_ID = "Plot_name",           
  draw_plot = TRUE, 
  rm_outliers = FALSE
)

st_write(
  obj = fieldplot_old$polygon,
  dsn = "~/Remote-sensing/Ranong_fieldplot_old.shp",
  delete_layer = TRUE
)

Deleting layer `Ranong_fieldplot_old' using driver `ESRI Shapefile'
Writing layer `Ranong_fieldplot_old' to data source 
  `/Users/pchwx/Remote-sensing/Ranong_fieldplot_old.shp' using driver `ESRI Shapefile'
Writing 1 features with 0 fields and geometry type Polygon.

st_write(
  obj = fieldplot_new$polygon,
  dsn = "~/Remote-sensing/Ranong_fieldplot_new.shp",
  delete_layer = TRUE
)

Deleting layer `Ranong_fieldplot_new' using driver `ESRI Shapefile'
Writing layer `Ranong_fieldplot_new' to data source 
  `/Users/pchwx/Remote-sensing/Ranong_fieldplot_new.shp' using driver `ESRI Shapefile'
Writing 1 features with 0 fields and geometry type Polygon.

# 20x20
subplots_grid20_old <- divide_plot(
  corner_data = fieldplot_old$corner_coord, 
  rel_coord = c("x_rel","y_rel"),
  proj_coord = c("x_proj","y_proj"),
  grid_size = 20,
  tree_data = data_forAGB, 
  tree_coords = c("PX","PY")
)

#50x50
subplots_grid50_old <- divide_plot(
  corner_data = fieldplot_old$corner_coord, 
  rel_coord = c("x_rel","y_rel"),
  proj_coord = c("x_proj","y_proj"),
  grid_size = 50,
  tree_data = data_forAGB, 
  tree_coords = c("PX","PY")
)

# 20x20
subplots_grid20_new <- divide_plot(
  corner_data = fieldplot_new$corner_coord, 
  rel_coord = c("x_rel","y_rel"),
  proj_coord = c("x_proj","y_proj"),
  grid_size = 20,
  tree_data = new_data_for_AGB, 
  tree_coords = c("PX","PY")
)

#50x50
subplots_grid50_new <- divide_plot(
  corner_data = fieldplot_old$corner_coord, 
  rel_coord = c("x_rel","y_rel"),
  proj_coord = c("x_proj","y_proj"),
  grid_size = 50,
  tree_data = data_forAGB, 
  tree_coords = c("PX","PY")
)

subplot_metric_perHA20_old <- subplot_summary(
  subplots = subplots_grid20_old,
  value = "AGB", 
  per_ha = TRUE)

[[1]]

subplot_metric_perHA20_new <- subplot_summary(
  subplots = subplots_grid20_new,
  value = "AGB", 
  per_ha = TRUE)

[[1]]

subplot_metric_perHA50_old <- subplot_summary(
  subplots = subplots_grid50_old,
  value = "AGB", 
  per_ha = TRUE)

[[1]]

subplot_metric_perHA50_new <- subplot_summary(
  subplots = subplots_grid50_new,
  value = "AGB", 
  per_ha = TRUE)

[[1]]

custom_plot_perHA20_old <- subplot_metric_perHA20_old$plot_design +
  scale_fill_distiller(palette = "YlGnBu", direction = 1) 

custom_plot_perHA20_old

sf::st_write(custom_plot_perHA20_old,"~/Remote-sensing/polygon_20m_perHA.shp")

Error in UseMethod("st_write") : 
  no applicable method for 'st_write' applied to an object of class "c('gg', 'ggplot')"

custom_plot_perHA20_new <- subplot_metric_perHA20_new$plot_design +
  scale_fill_distiller(palette = "YlGnBu", direction = 1)

custom_plot_perHA20_new

subplot_polygons20m_perHA <- sf::st_set_crs(
 subplot_metric_perHA20_new$polygon ,
  value = "EPSG:32647") 

sf::st_write(subplot_polygons20m_perHA,"~/Remote-sensing/polygon_20m_perHA.shp")

Writing layer `polygon_20m_perHA' to data source `/Users/pchwx/Remote-sensing/polygon_20m_perHA.shp' using driver `ESRI Shapefile'
Writing 225 features with 2 fields and geometry type Polygon.

custom_plot_perHA50_old <- subplot_metric_perHA50_old$plot_design +
  scale_fill_distiller(palette = "YlGnBu", direction = 1) 

custom_plot_perHA50_old

custom_plot_perHA50_new <- subplot_metric_perHA50_new$plot_design +
  scale_fill_distiller(palette = "YlGnBu", direction = 1)

custom_plot_perHA50_new

subplot_metric_sum20_old<- subplot_summary(
  subplots = subplots_grid20_old,
  value = "AGB",
  fun = quantile, probs = 0.5, 
  per_ha = FALSE)

[[1]]

custom_plot_qt <- subplot_metric_sum20_old$plot_design +
  scale_fill_distiller(palette = "YlGnBu", direction = 1)

custom_plot_qt